### Model fit exercise: calculate the estimated probability of war for the observations in our
### sample that actually went to war

date()
library("tidyverse")
library("foreach")
source("backend_main.r")
sessionInfo()


load("results/fit_base.rda")


## Calculate average (across imputed models/data) estimated probability of war for each dispute in
## the data
pr_war <- foreach (imp = fit_base, .combine = "cbind") %do% {
  setup <- structwar_setup(f_dispute = imp$f_dispute,
                           f_participant = imp$f_participant,
                           data_dispute = imp$data_dispute,
                           data_participant = imp$data_participant,
                           xlev_dispute = imp$xlev_dispute,
                           xlev_participant = imp$xlev_participant,
                           n_halton = imp$n_halton)
  parms <- extract_params(est = coef(imp$fit),
                          setup = setup,
                          for_counterfactuals = TRUE)
  cvals <- contest_vals(dispute_level_id = setup$dispute_level_id,
                        state_level_id = setup$state_level_id,
                        side_a = setup$side_a,
                        ratio = parms$ratio,
                        xmax = .Machine$double.xmax)
  prob_war(p_theta_a = setup$p_theta_a,
           loc_a = parms$loc_a,
           loc_b = parms$loc_b,
           scl_a = parms$scl_a,
           scl_b = parms$scl_b,
           pi_a = cvals$pi_a,
           pi_b = cvals$pi_b)
}
pr_war <- rowMeans(pr_war)
stopifnot(length(pr_war) == nrow(fit_base[[1]]$data_dispute))

## Extract probabilities associated with cases that actually went to war
data_wars <- fit_base[[1]]$data_dispute %>%
  add_column(pr_war_est = pr_war) %>%
  filter(war == 1) %>%
  select(id, n_states_a, n_states_b, pr_war_est) %>%
  arrange(desc(pr_war_est))

## Nice names for top and bottom wars
data_war_names <- tribble(
  ~id, ~war_name,
  3957, "Gulf War",
  51, "Korean War",
  257, "World War I",
  611, "Vietnam War",
  258, "World War II",
  349, "Sino-Soviet Border Conflict",
  189, "Russo-Turkish War",
  2749, "Sino-Vietnamese Border Conflict",
  3628, "Sino-Vietnamese Border Conflict",
  3622, "Sino-Vietnamese Border Conflict",
  3872, "Operation München",
  1163, "Acre War",
  1206, "Football War",
  1646, "French Invasion of Amapá",
  1279, "Bloody Christmas",
  1773, "Portugual/Germany in Angola"
)

## Clean up participant data to help make sense of which wars these are
data_war_info <- fit_base[[1]]$data_participant %>%
  group_by(id) %>%
  summarise(year = min(year),
            a_max = stateabb[sidea == 1 & cinc == max(cinc[sidea == 1])],
            b_max = stateabb[sidea == 0 & cinc == max(cinc[sidea == 0])])
data_wars <- data_wars %>%
  left_join(data_war_info, by = "id") %>%
  left_join(data_war_names, by = "id") %>%
  mutate(rank = row_number()) %>%
  select(rank, id, war_name, year, pr_war_est)

## Separate the tables to save space in JOP submission
tab_out_top <- c(
  "\\begin{tabular}{lllc}",
  "\\toprule",
  "Rank & Description & Year & $\\Pr$(War) \\\\",
  "\\midrule",
  with(data_wars %>% head(5),
       str_glue("{rank} & {war_name} & {year} & {sprintf('%.2f', pr_war_est)} \\\\")),
  "\\bottomrule",
  "\\end{tabular}"
)
tab_out_bottom <- c(
  "\\begin{tabular}{lllc}",
  "\\toprule",
  "Rank & Description & Year & $\\Pr$(War) \\\\",
  "\\midrule",
  with(data_wars %>% tail(5),
       str_glue("{rank} & {war_name} & {year} & {sprintf('%.2f', pr_war_est)} \\\\")),
  "\\bottomrule",
  "\\end{tabular}"
)

if (!dir.exists("tables"))
  dir.create("tables")
writeLines(tab_out_top, con = "tables/table_4a.tex")
writeLines(tab_out_bottom, con = "tables/table_4b.tex")


date()
